---
title: "Relational Support and School Belonging as Indicators of Childrens Subjective WellBeing in the UK"
author: "Tendai Charles"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
---

## Project overview

This R Markdown file estimates and compares SEMs for the MiCREATE UK survey.

**Research question:**  
To what extent do family support and peer support contribute to children's subjective wellbeing directly and indirectly through school belonging?

### What this R script does

This script estimates:

1. A **CFA** for subjective wellbeing and school belonging  
2. A **baseline partial-mediation SEM**  
3. A **full-mediation SEM**  
4. A **reliability-corrected SEM**  
5. A **refined reliability-corrected SEM** with `Q24a ~~ Q24b`

It also produces:

- a **clean structural SEM path diagram** 
- a **full measurement-model diagram** 

### Practical note

This file assumes that the `.sav` file is in the **same folder** as this `.Rmd` file.

The expected filename is:

`MiCREATE_WP567_UK_MiCREATE_Survey_data.sav`

If you are replicating my study, but your file has a different name, edit the `data_path` object in Section 2.

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)

options(width = 120, scipen = 999)
set.seed(20260312)
```

## 1. Install and load packages

```{r packages}
required_packages <- c(
  "haven",
  "dplyr",
  "tidyr",
  "tibble",
  "knitr",
  "psych",
  "lavaan",
  "semPlot"
)

to_install <- setdiff(required_packages, rownames(installed.packages()))
if (length(to_install) > 0) {
  install.packages(to_install, repos = "https://cloud.r-project.org", dependencies = TRUE)
}

invisible(
  lapply(required_packages, function(pkg) {
    library(pkg, character.only = TRUE)
  })
)
```

## 2. Read the data

```{r load-data}
data_path <- "MiCREATE_WP567_UK_MiCREATE_Survey_data.sav"

if (!file.exists(data_path)) {
  stop(
    paste0(
      "The SPSS file was not found at: ", data_path, "\n",
      "Please either:\n",
      "1. place the .sav file in the same folder as this .Rmd file, or\n",
      "2. edit the 'data_path' object in Section 2."
    )
  )
}

data_raw <- haven::read_sav(data_path, user_na = FALSE)

dim(data_raw)
names(data_raw)[1:min(20, ncol(data_raw))]
```

## 3. Define project variables

```{r define-variables}
wellbeing_items <- c("Q24a", "Q24b", "Q24c", "Q24d")
family_items    <- c("Q29a", "Q29b")
peer_items      <- c("Q29c", "Q29d")
belonging_items <- c("Q34a", "Q34b", "Q34c", "Q34d", "Q34e")

all_project_items <- c(wellbeing_items, family_items, peer_items, belonging_items)

missing_vars <- setdiff(all_project_items, names(data_raw))

if (length(missing_vars) > 0) {
  stop(
    paste(
      "The following expected variables were not found in the .sav file:",
      paste(missing_vars, collapse = ", ")
    )
  )
}

all_project_items
```

## 4. Helper functions

```{r helper-functions}
lavaan_estimator <- "MLR"
missing_handling <- "fiml"

get_label <- function(x) {
  lab <- attr(x, "label")
  if (is.null(lab)) NA_character_ else as.character(lab)
}

recode_project_item <- function(x) {
  x <- as.numeric(x)
  x[x < 0] <- NA_real_
  x[x < 1 | x > 5] <- NA_real_
  x
}

scale_mean <- function(dat, vars, min_valid = 1) {
  x <- dat[, vars, drop = FALSE]
  out <- rowMeans(x, na.rm = TRUE)
  out[rowSums(!is.na(x)) < min_valid] <- NA_real_
  out
}

safe_alpha <- function(dat) {
  out <- tryCatch(
    psych::alpha(dat),
    error = function(e) e
  )
  if (inherits(out, "error")) NA_real_ else out$total$raw_alpha
}

safe_variance <- function(x) {
  if (sum(!is.na(x)) <= 1) return(NA_real_)
  stats::var(x, na.rm = TRUE)
}

safe_error_variance <- function(var_obs, alpha_value, fallback_alpha = 0.70) {
  if (is.na(alpha_value) || alpha_value <= 0 || alpha_value >= 1) {
    alpha_value <- fallback_alpha
  }
  out <- var_obs * (1 - alpha_value)
  if (is.na(out) || out <= 0) out <- var_obs * 0.20
  out
}

is_converged <- function(fit) {
  inherits(fit, "lavaan") &&
    isTRUE(
      tryCatch(
        lavaan::lavInspect(fit, "converged"),
        error = function(e) FALSE
      )
    )
}

model_note <- function(fit) {
  if (inherits(fit, "error")) {
    return(paste("Estimation error:", conditionMessage(fit)))
  }
  if (!inherits(fit, "lavaan")) {
    return("Object is not a lavaan model.")
  }
  if (!is_converged(fit)) {
    return("Model did not converge.")
  }
  NA_character_
}

safe_fit_value <- function(fit, candidates) {
  if (!is_converged(fit)) return(NA_real_)
  fm <- lavaan::fitMeasures(fit)
  hit <- candidates[candidates %in% names(fm)][1]
  if (length(hit) == 0 || is.na(hit)) return(NA_real_)
  unname(fm[hit])
}

extract_fit_table <- function(fit) {
  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        fit_index = c(
          "chisq.scaled", "df.scaled", "pvalue.scaled",
          "cfi.robust", "tli.robust", "rmsea.robust",
          "srmr", "aic", "bic"
        ),
        value = NA_real_,
        note = model_note(fit)
      )
    )
  }

  fm <- lavaan::fitMeasures(fit)

  wanted <- c(
    if ("chisq.scaled" %in% names(fm)) "chisq.scaled" else "chisq",
    if ("df.scaled" %in% names(fm)) "df.scaled" else "df",
    if ("pvalue.scaled" %in% names(fm)) "pvalue.scaled" else "pvalue",
    if ("cfi.robust" %in% names(fm)) "cfi.robust" else "cfi",
    if ("tli.robust" %in% names(fm)) "tli.robust" else "tli",
    if ("rmsea.robust" %in% names(fm)) "rmsea.robust" else "rmsea",
    if ("srmr" %in% names(fm)) "srmr" else NULL,
    if ("aic" %in% names(fm)) "aic" else NULL,
    if ("bic" %in% names(fm)) "bic" else NULL
  )

  wanted <- unique(wanted)

  tibble::tibble(
    fit_index = wanted,
    value = unname(lavaan::fitMeasures(fit, fit.measures = wanted)),
    note = NA_character_
  )
}

extract_parameter_table <- function(fit, ops = c("=~", "~", "~~", ":=")) {
  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        lhs = character(),
        op = character(),
        rhs = character(),
        est = numeric(),
        se = numeric(),
        z = numeric(),
        pvalue = numeric(),
        ci.lower = numeric(),
        ci.upper = numeric(),
        std.all = numeric()
      )
    )
  }

  lavaan::parameterEstimates(
    fit,
    standardized = TRUE,
    ci = TRUE
  ) %>%
    dplyr::filter(op %in% ops) %>%
    dplyr::select(lhs, op, rhs, est, se, z, pvalue, ci.lower, ci.upper, std.all)
}

extract_key_paths_baseline <- function(fit) {
  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        lhs = character(),
        op = character(),
        rhs = character(),
        est = numeric(),
        se = numeric(),
        z = numeric(),
        pvalue = numeric(),
        ci.lower = numeric(),
        ci.upper = numeric(),
        std.all = numeric()
      )
    )
  }

  pe <- lavaan::parameterEstimates(
    fit,
    standardized = TRUE,
    ci = TRUE
  )

  direct_paths <- pe %>%
    dplyr::filter(
      op == "~",
      (
        (lhs == "school_belonging" & rhs %in% c("family_support", "peer_support")) |
          (lhs == "wellbeing" & rhs %in% c("school_belonging", "family_support", "peer_support"))
      )
    )

  indirect_paths <- pe %>%
    dplyr::filter(
      op == ":=",
      lhs %in% c("ind_family", "ind_peer", "total_family", "total_peer")
    )

  dplyr::bind_rows(direct_paths, indirect_paths) %>%
    dplyr::select(lhs, op, rhs, est, se, z, pvalue, ci.lower, ci.upper, std.all)
}

extract_key_paths_rc <- function(fit) {
  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        lhs = character(),
        op = character(),
        rhs = character(),
        est = numeric(),
        se = numeric(),
        z = numeric(),
        pvalue = numeric(),
        ci.lower = numeric(),
        ci.upper = numeric(),
        std.all = numeric()
      )
    )
  }

  pe <- lavaan::parameterEstimates(
    fit,
    standardized = TRUE,
    ci = TRUE
  )

  direct_paths <- pe %>%
    dplyr::filter(
      op == "~",
      (
        (lhs == "school_belonging" & rhs %in% c("family_lat", "peer_lat")) |
          (lhs == "wellbeing" & rhs %in% c("school_belonging", "family_lat", "peer_lat"))
      )
    )

  indirect_paths <- pe %>%
    dplyr::filter(
      op == ":=",
      lhs %in% c("ind_family", "ind_peer", "total_family", "total_peer")
    )

  dplyr::bind_rows(direct_paths, indirect_paths) %>%
    dplyr::select(lhs, op, rhs, est, se, z, pvalue, ci.lower, ci.upper, std.all)
}

extract_model_row <- function(fit, label) {
  tibble::tibble(
    model = label,
    chisq_scaled = safe_fit_value(fit, c("chisq.scaled", "chisq")),
    df_scaled = safe_fit_value(fit, c("df.scaled", "df")),
    pvalue_scaled = safe_fit_value(fit, c("pvalue.scaled", "pvalue")),
    cfi_robust = safe_fit_value(fit, c("cfi.robust", "cfi")),
    tli_robust = safe_fit_value(fit, c("tli.robust", "tli")),
    rmsea_robust = safe_fit_value(fit, c("rmsea.robust", "rmsea")),
    srmr = safe_fit_value(fit, c("srmr")),
    aic = safe_fit_value(fit, c("aic")),
    bic = safe_fit_value(fit, c("bic")),
    note = model_note(fit)
  )
}

extract_nested_test <- function(fit_restricted,
                                fit_full,
                                restricted_label = "restricted",
                                full_label = "full",
                                method = "satorra.bentler.2001") {
  if (!is_converged(fit_restricted) || !is_converged(fit_full)) {
    return(
      tibble::tibble(
        restricted_model = restricted_label,
        full_model = full_label,
        chisq_diff = NA_real_,
        df_diff = NA_real_,
        p_value = NA_real_,
        note = "At least one model did not converge."
      )
    )
  }

  lrt <- tryCatch(
    lavaan::lavTestLRT(fit_restricted, fit_full, method = method),
    error = function(e) e
  )

  if (inherits(lrt, "error")) {
    return(
      tibble::tibble(
        restricted_model = restricted_label,
        full_model = full_label,
        chisq_diff = NA_real_,
        df_diff = NA_real_,
        p_value = NA_real_,
        note = conditionMessage(lrt)
      )
    )
  }

  lrt_df <- as.data.frame(lrt)

  chisq_diff_col <- names(lrt_df)[grepl("Chisq diff", names(lrt_df), fixed = TRUE)][1]
  df_diff_col    <- names(lrt_df)[grepl("Df diff", names(lrt_df), fixed = TRUE)][1]
  p_col          <- names(lrt_df)[grepl("^Pr\\(>Chisq\\)$", names(lrt_df))][1]

  candidate_rows <- integer(0)

  if (!is.na(chisq_diff_col) && length(chisq_diff_col) == 1) {
    candidate_rows <- which(!is.na(lrt_df[[chisq_diff_col]]))
  }

  if (length(candidate_rows) == 0 && !is.na(p_col) && length(p_col) == 1) {
    candidate_rows <- which(!is.na(lrt_df[[p_col]]))
  }

  if (length(candidate_rows) == 0) {
    candidate_rows <- nrow(lrt_df)
  }

  i <- candidate_rows[length(candidate_rows)]

  tibble::tibble(
    restricted_model = restricted_label,
    full_model = full_label,
    chisq_diff = if (!is.na(chisq_diff_col) && length(chisq_diff_col) == 1) as.numeric(lrt_df[[chisq_diff_col]][i]) else NA_real_,
    df_diff = if (!is.na(df_diff_col) && length(df_diff_col) == 1) as.numeric(lrt_df[[df_diff_col]][i]) else NA_real_,
    p_value = if (!is.na(p_col) && length(p_col) == 1) as.numeric(lrt_df[[p_col]][i]) else NA_real_,
    note = NA_character_
  )
}

get_std_est <- function(fit, lhs, op, rhs) {
  if (!is_converged(fit)) return(NA_real_)

  pe <- lavaan::parameterEstimates(fit, standardized = TRUE)
  out <- pe$std.all[pe$lhs == lhs & pe$op == op & pe$rhs == rhs]

  if (length(out) == 0) return(NA_real_)
  as.numeric(out[1])
}

get_full_semplot_labels <- function(fit) {
  tryCatch({
    spm <- semPlot::semPlotModel(fit)
    labs <- spm@Vars$name

    pretty_map <- c(
      family_lat = "Family\nsupport",
      peer_lat = "Peer\nsupport",
      school_belonging = "School\nbelonging",
      wellbeing = "Subjective\nwellbeing",
      family_support = "Family\nsupport\n(composite)",
      peer_support = "Peer\nsupport\n(composite)",
      Q24a = "Satisfied\nwith life",
      Q24b = "Have what I\nwant in life",
      Q24c = "Future looks\ngood",
      Q24d = "Feel good\nabout life",
      Q34a = "Like being\nin school",
      Q34b = "Feel safe\nat school",
      Q34c = "Belong in\nthis class",
      Q34d = "Comfortable\nasking teachers",
      Q34e = "Okay with\nclassmates"
    )

    idx <- labs %in% names(pretty_map)
    labs[idx] <- pretty_map[labs[idx]]
    labs
  }, error = function(e) NULL)
}

save_full_measurement_diagram <- function(fit,
                                          file_stub,
                                          diagram_dir,
                                          width_px = 4200,
                                          height_px = 3200,
                                          width_in = 14,
                                          height_in = 10) {
  png_file <- file.path(diagram_dir, paste0(file_stub, ".png"))
  pdf_file <- file.path(diagram_dir, paste0(file_stub, ".pdf"))

  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        file = c(png_file, pdf_file),
        created = FALSE,
        note = "Model did not converge; full measurement diagram not created."
      )
    )
  }

  node_labels <- get_full_semplot_labels(fit)
  draw_mode <- NA_character_

  draw_semplot <- function(with_labels = TRUE) {
    args <- list(
      object = fit,
      what = "std",
      whatLabels = "std",
      style = "lisrel",
      layout = "tree",
      rotation = 2,
      intercepts = FALSE,
      residuals = FALSE,
      exoCov = TRUE,
      curvePivot = TRUE,
      edge.label.cex = 0.90,
      sizeLat = 11,
      sizeMan = 7,
      nCharNodes = 0,
      mar = c(8, 8, 8, 8),
      title = FALSE
    )

    if (with_labels && !is.null(node_labels)) {
      args$nodeLabels <- node_labels
    }

    do.call(semPlot::semPaths, args)
  }

  render_plot <- function() {
    tryCatch({
      draw_semplot(with_labels = TRUE)
      draw_mode <<- "Created with readable node labels."
    }, error = function(e) {
      draw_semplot(with_labels = FALSE)
      draw_mode <<- "Created with default node labels."
    })
  }

  png_ok <- tryCatch({
    grDevices::png(
      filename = png_file,
      width = width_px,
      height = height_px,
      res = 300,
      bg = "white"
    )
    render_plot()
    grDevices::dev.off()
    TRUE
  }, error = function(e) {
    try(grDevices::dev.off(), silent = TRUE)
    FALSE
  })

  pdf_ok <- tryCatch({
    grDevices::pdf(
      file = pdf_file,
      width = width_in,
      height = height_in,
      useDingbats = FALSE
    )
    render_plot()
    grDevices::dev.off()
    TRUE
  }, error = function(e) {
    try(grDevices::dev.off(), silent = TRUE)
    FALSE
  })

  tibble::tibble(
    file = c(png_file, pdf_file),
    created = c(png_ok, pdf_ok),
    note = c(draw_mode, draw_mode)
  )
}


save_structural_path_diagram <- function(fit,
                                         file_stub,
                                         diagram_dir,
                                         width_px = 3800,
                                         height_px = 2200,
                                         width_in = 12.5,
                                         height_in = 7.2) {
  png_file <- file.path(diagram_dir, paste0(file_stub, ".png"))
  pdf_file <- file.path(diagram_dir, paste0(file_stub, ".pdf"))

  if (!is_converged(fit)) {
    return(
      tibble::tibble(
        file = c(png_file, pdf_file),
        created = FALSE,
        note = "Model did not converge; structural diagram not created."
      )
    )
  }

  b_family_belong <- get_std_est(fit, "school_belonging", "~", "family_lat")
  b_peer_belong   <- get_std_est(fit, "school_belonging", "~", "peer_lat")
  b_belong_well   <- get_std_est(fit, "wellbeing", "~", "school_belonging")
  b_family_well   <- get_std_est(fit, "wellbeing", "~", "family_lat")
  b_peer_well     <- get_std_est(fit, "wellbeing", "~", "peer_lat")
  cov_family_peer <- get_std_est(fit, "family_lat", "~~", "peer_lat")

  r2_vals <- tryCatch(
    lavaan::inspect(fit, "r2"),
    error = function(e) numeric(0)
  )

  r2_school <- if ("school_belonging" %in% names(r2_vals)) as.numeric(r2_vals["school_belonging"]) else NA_real_
  r2_well   <- if ("wellbeing" %in% names(r2_vals)) as.numeric(r2_vals["wellbeing"]) else NA_real_

  draw_ellipse <- function(x, y, a = 1.20, b = 0.65, lwd = 2) {
    theta <- seq(0, 2 * pi, length.out = 400)
    lines(x + a * cos(theta), y + b * sin(theta), lwd = lwd)
  }

  draw_beta_label <- function(x, y, val, cex = 0.96) {
    text(
      x, y,
      labels = bquote(beta == .(sprintf("%.3f", val))),
      cex = cex
    )
  }

  draw_r_label <- function(x, y, val, cex = 0.94) {
    text(
      x, y,
      labels = bquote(italic(r) == .(sprintf("%.3f", val))),
      cex = cex
    )
  }

  draw_r2_label <- function(x, y, val, cex = 0.92) {
    if (!is.na(val)) {
      text(
        x, y,
        labels = bquote(R^2 == .(sprintf("%.3f", val))),
        cex = cex
      )
    }
  }

  draw_curved_double_arrow <- function(start,
                                       control,
                                       end,
                                       lwd = 2,
                                       length = 0.09,
                                       n = 200) {
    t <- seq(0, 1, length.out = n)

    x <- (1 - t)^2 * start[1] + 2 * (1 - t) * t * control[1] + t^2 * end[1]
    y <- (1 - t)^2 * start[2] + 2 * (1 - t) * t * control[2] + t^2 * end[2]

    lines(x, y, lwd = lwd)

    arrows(
      x0 = x[2], y0 = y[2],
      x1 = x[1], y1 = y[1],
      length = length,
      lwd = lwd
    )

    arrows(
      x0 = x[n - 1], y0 = y[n - 1],
      x1 = x[n], y1 = y[n],
      length = length,
      lwd = lwd
    )
  }

  draw_structural <- function() {
    old_par <- par(no.readonly = TRUE)
    on.exit(par(old_par), add = TRUE)

    par(mar = c(0.35, 0.35, 0.35, 0.35), xpd = NA)

    plot(
      NA, NA,
      type = "n",
      xlim = c(0.35, 10.15),
      ylim = c(0.95, 5.05),
      axes = FALSE,
      xlab = "",
      ylab = "",
      bty = "n"
    )

    family_xy <- c(2.15, 4.00)
    peer_xy   <- c(2.15, 1.95)
    sb_xy     <- c(5.75, 2.98)
    wb_xy     <- c(9.10, 2.98)

    ex_a  <- 1.10
    ex_b  <- 0.56
    lat_a <- 1.32
    lat_b <- 0.70

    draw_ellipse(family_xy[1], family_xy[2], ex_a, ex_b, lwd = 2)
    draw_ellipse(peer_xy[1], peer_xy[2], ex_a, ex_b, lwd = 2)
    draw_ellipse(sb_xy[1], sb_xy[2], lat_a, lat_b, lwd = 2)
    draw_ellipse(wb_xy[1], wb_xy[2], lat_a, lat_b, lwd = 2)

    text(family_xy[1], family_xy[2], "Family\nsupport", cex = 1.10, font = 2)
    text(peer_xy[1], peer_xy[2], "Peer\nsupport", cex = 1.10, font = 2)

    text(sb_xy[1], sb_xy[2] + 0.14, "School\nbelonging", cex = 1.06, font = 2)
    draw_r2_label(sb_xy[1], sb_xy[2] - 0.46, r2_school, cex = 0.92)

    text(wb_xy[1], wb_xy[2] + 0.14, "Subjective\nwellbeing", cex = 1.06, font = 2)
    draw_r2_label(wb_xy[1], wb_xy[2] - 0.46, r2_well, cex = 0.92)

    draw_curved_double_arrow(
      start = c(family_xy[1] - ex_a * 0.86, family_xy[2] - 0.14),
      control = c(0.75, 2.98),
      end = c(peer_xy[1] - ex_a * 0.86, peer_xy[2] + 0.14),
      lwd = 1.9,
      length = 0.08
    )
    draw_r_label(0.95, 2.98, cov_family_peer, cex = 0.92)

    arrows(
      x0 = family_xy[1] + ex_a * 0.92, y0 = family_xy[2] - 0.05,
      x1 = sb_xy[1] - lat_a * 0.95,   y1 = sb_xy[2] + 0.24,
      length = 0.085,
      lwd = 2
    )
    draw_beta_label(4.00, 3.72, b_family_belong, cex = 0.94)

    arrows(
      x0 = peer_xy[1] + ex_a * 0.92, y0 = peer_xy[2] + 0.05,
      x1 = sb_xy[1] - lat_a * 0.95,  y1 = sb_xy[2] - 0.24,
      length = 0.085,
      lwd = 2
    )
    draw_beta_label(4.00, 2.22, b_peer_belong, cex = 0.94)

    arrows(
      x0 = sb_xy[1] + lat_a * 0.98, y0 = sb_xy[2],
      x1 = wb_xy[1] - lat_a * 0.98, y1 = wb_xy[2],
      length = 0.085,
      lwd = 2
    )
    draw_beta_label(7.42, 3.33, b_belong_well, cex = 0.94)

    arrows(
      x0 = family_xy[1] + ex_a * 0.92, y0 = family_xy[2] + 0.18,
      x1 = wb_xy[1] - lat_a * 0.96,    y1 = wb_xy[2] + 0.34,
      length = 0.085,
      lwd = 2
    )
    draw_beta_label(5.65, 4.36, b_family_well, cex = 0.94)

    arrows(
      x0 = peer_xy[1] + ex_a * 0.92, y0 = peer_xy[2] - 0.18,
      x1 = wb_xy[1] - lat_a * 0.96,  y1 = wb_xy[2] - 0.34,
      length = 0.085,
      lwd = 2
    )
    draw_beta_label(5.72, 1.76, b_peer_well, cex = 0.94)
  }

  png_ok <- tryCatch({
    grDevices::png(
      filename = png_file,
      width = width_px,
      height = height_px,
      res = 300,
      bg = "white"
    )
    draw_structural()
    grDevices::dev.off()
    TRUE
  }, error = function(e) {
    try(grDevices::dev.off(), silent = TRUE)
    FALSE
  })

  pdf_ok <- tryCatch({
    grDevices::pdf(
      file = pdf_file,
      width = width_in,
      height = height_in,
      useDingbats = FALSE
    )
    draw_structural()
    grDevices::dev.off()
    TRUE
  }, error = function(e) {
    try(grDevices::dev.off(), silent = TRUE)
    FALSE
  })

  tibble::tibble(
    file = c(png_file, pdf_file),
    created = c(png_ok, pdf_ok),
    note = c(
      "Journal-style structural diagram with ellipses, curved covariance arrow, tighter spacing, and R-squared values.",
      "Journal-style structural diagram with ellipses, curved covariance arrow, tighter spacing, and R-squared values."
    )
  )
}

```

## 5. Prepare the analytic dataset

```{r prepare-data}
data_sem <- data_raw %>%
  dplyr::select(dplyr::all_of(all_project_items)) %>%
  dplyr::mutate(
    dplyr::across(dplyr::all_of(all_project_items), recode_project_item)
  )

data_sem$family_support <- scale_mean(data_sem, family_items, min_valid = 1)
data_sem$peer_support   <- scale_mean(data_sem, peer_items, min_valid = 1)

data_sem$wellbeing_mean        <- scale_mean(data_sem, wellbeing_items, min_valid = 2)
data_sem$school_belonging_mean <- scale_mean(data_sem, belonging_items, min_valid = 3)

analysis_variables <- c(
  all_project_items,
  "family_support",
  "peer_support",
  "wellbeing_mean",
  "school_belonging_mean"
)

str(data_sem[, analysis_variables, drop = FALSE])
```

## 6. Item dictionary

```{r item-dictionary}
item_dictionary <- tibble::tibble(
  variable = all_project_items,
  label = sapply(data_raw[all_project_items], get_label)
)

knitr::kable(item_dictionary, caption = "Item dictionary for the project")
```

## 7. Descriptives, missingness, and reliability

```{r descriptives}
descriptives <- psych::describe(data_sem[, analysis_variables, drop = FALSE]) %>%
  as.data.frame() %>%
  tibble::rownames_to_column("variable")

knitr::kable(
  descriptives[, c("variable", "n", "mean", "sd", "min", "max", "skew", "kurtosis")],
  digits = 2,
  caption = "Descriptive statistics"
)
```

```{r missingness}
missingness <- data_sem %>%
  dplyr::summarise(
    dplyr::across(
      dplyr::all_of(analysis_variables),
      ~ sum(is.na(.))
    )
  ) %>%
  tidyr::pivot_longer(
    cols = dplyr::everything(),
    names_to = "variable",
    values_to = "n_missing"
  ) %>%
  dplyr::mutate(
    pct_missing = 100 * n_missing / nrow(data_sem)
  )

knitr::kable(missingness, digits = 2, caption = "Missingness by variable")
```

```{r reliability}
alpha_wellbeing <- safe_alpha(data_sem[, wellbeing_items])
alpha_family    <- safe_alpha(data_sem[, family_items])
alpha_peer      <- safe_alpha(data_sem[, peer_items])
alpha_belonging <- safe_alpha(data_sem[, belonging_items])

reliability_table <- tibble::tibble(
  scale = c("Subjective wellbeing", "Family support", "Peer support", "School belonging"),
  alpha = c(alpha_wellbeing, alpha_family, alpha_peer, alpha_belonging)
)

knitr::kable(reliability_table, digits = 3, caption = "Cronbach's alpha for project scales")
```

## 8. Measurement model (CFA)

```{r cfa-model}
measurement_model <- '
  wellbeing =~ Q24a + Q24b + Q24c + Q24d
  school_belonging =~ Q34a + Q34b + Q34c + Q34d + Q34e
'

fit_cfa <- tryCatch(
  lavaan::cfa(
    model = measurement_model,
    data = data_sem,
    estimator = lavaan_estimator,
    std.lv = TRUE,
    missing = missing_handling
  ),
  error = function(e) e
)
```

```{r cfa-fit}
cfa_fit_table <- extract_fit_table(fit_cfa)
knitr::kable(cfa_fit_table, digits = 4, caption = "CFA fit indices")
```

```{r cfa-loadings}
cfa_loadings <- extract_parameter_table(fit_cfa) %>%
  dplyr::filter(op == "=~")

if (nrow(cfa_loadings) > 0) {
  knitr::kable(cfa_loadings, digits = 3, caption = "CFA standardized loadings")
} else {
  cat("No CFA loadings are available because the CFA did not converge.\n")
}
```

## 9. Baseline structural models

### 9.1 Partial mediation model

```{r sem-partial-model}
sem_partial_model <- '
  wellbeing =~ Q24a + Q24b + Q24c + Q24d
  school_belonging =~ Q34a + Q34b + Q34c + Q34d + Q34e

  school_belonging ~ a1*family_support + a2*peer_support
  wellbeing ~ b*school_belonging + c1*family_support + c2*peer_support

  family_support ~~ peer_support

  ind_family := a1*b
  ind_peer   := a2*b
  total_family := c1 + (a1*b)
  total_peer   := c2 + (a2*b)
'

fit_sem_partial <- tryCatch(
  lavaan::sem(
    model = sem_partial_model,
    data = data_sem,
    estimator = lavaan_estimator,
    std.lv = TRUE,
    fixed.x = FALSE,
    missing = missing_handling
  ),
  error = function(e) e
)
```

```{r sem-partial-fit}
sem_partial_fit <- extract_fit_table(fit_sem_partial)
knitr::kable(sem_partial_fit, digits = 4, caption = "Baseline partial-mediation SEM fit indices")
```

```{r sem-partial-key}
sem_partial_key <- extract_key_paths_baseline(fit_sem_partial)

if (nrow(sem_partial_key) > 0) {
  knitr::kable(sem_partial_key, digits = 3, caption = "Baseline partial-mediation SEM: key structural paths and effects")
} else {
  cat("No baseline partial-mediation parameter table is available because the model did not converge.\n")
}
```

### 9.2 Full mediation model

```{r sem-full-model}
sem_full_model <- '
  wellbeing =~ Q24a + Q24b + Q24c + Q24d
  school_belonging =~ Q34a + Q34b + Q34c + Q34d + Q34e

  school_belonging ~ a1*family_support + a2*peer_support
  wellbeing ~ b*school_belonging

  family_support ~~ peer_support

  ind_family := a1*b
  ind_peer   := a2*b
'

fit_sem_full <- tryCatch(
  lavaan::sem(
    model = sem_full_model,
    data = data_sem,
    estimator = lavaan_estimator,
    std.lv = TRUE,
    fixed.x = FALSE,
    missing = missing_handling
  ),
  error = function(e) e
)
```

```{r sem-full-fit}
sem_full_fit <- extract_fit_table(fit_sem_full)
knitr::kable(sem_full_fit, digits = 4, caption = "Full-mediation SEM fit indices")
```

## 10. Corrected comparison: full mediation versus partial mediation

```{r full-vs-partial-model-summary}
full_vs_partial_models <- dplyr::bind_rows(
  extract_model_row(fit_sem_full, "Full mediation"),
  extract_model_row(fit_sem_partial, "Partial mediation")
)

knitr::kable(full_vs_partial_models, digits = 4, caption = "Model-wise fit summary: full mediation and partial mediation")
```

```{r full-vs-partial-diff-test}
full_vs_partial_test <- extract_nested_test(
  fit_restricted = fit_sem_full,
  fit_full = fit_sem_partial,
  restricted_label = "Full mediation",
  full_label = "Partial mediation"
)

knitr::kable(full_vs_partial_test, digits = 4, caption = "Difference test: full mediation versus partial mediation")
```

## 11. Reliability-correction constants

```{r reliability-correction}
var_family_support <- safe_variance(data_sem$family_support)
var_peer_support   <- safe_variance(data_sem$peer_support)

err_var_family <- safe_error_variance(var_family_support, alpha_family)
err_var_peer   <- safe_error_variance(var_peer_support, alpha_peer)

reliability_correction_table <- tibble::tibble(
  composite = c("family_support", "peer_support"),
  observed_variance = c(var_family_support, var_peer_support),
  alpha_used = c(alpha_family, alpha_peer),
  fixed_error_variance = c(err_var_family, err_var_peer)
)

knitr::kable(reliability_correction_table, digits = 4, caption = "Reliability-correction constants used in the single-indicator latent models")
```

## 12. Reliability-corrected SEM (base model)

```{r rc-base-model}
sem_rc_base_lines <- c(
  "family_lat =~ 1*family_support",
  "peer_lat =~ 1*peer_support",
  paste0("family_support ~~ ", format(err_var_family, digits = 8, nsmall = 6), "*family_support"),
  paste0("peer_support ~~ ", format(err_var_peer, digits = 8, nsmall = 6), "*peer_support"),
  "",
  "wellbeing =~ Q24a + Q24b + Q24c + Q24d",
  "school_belonging =~ Q34a + Q34b + Q34c + Q34d + Q34e",
  "",
  "school_belonging ~ a1*family_lat + a2*peer_lat",
  "wellbeing ~ b*school_belonging + c1*family_lat + c2*peer_lat",
  "",
  "family_lat ~~ peer_lat",
  "",
  "ind_family := a1*b",
  "ind_peer   := a2*b",
  "total_family := c1 + (a1*b)",
  "total_peer   := c2 + (a2*b)"
)

sem_rc_base_model <- paste(sem_rc_base_lines, collapse = "\n")

fit_sem_rc_base <- tryCatch(
  lavaan::sem(
    model = sem_rc_base_model,
    data = data_sem,
    estimator = lavaan_estimator,
    fixed.x = FALSE,
    missing = missing_handling
  ),
  error = function(e) e
)
```

```{r rc-base-fit}
sem_rc_base_fit <- extract_fit_table(fit_sem_rc_base)
knitr::kable(sem_rc_base_fit, digits = 4, caption = "Reliability-corrected SEM fit indices")
```

```{r rc-base-key}
sem_rc_base_key <- extract_key_paths_rc(fit_sem_rc_base)

if (nrow(sem_rc_base_key) > 0) {
  knitr::kable(sem_rc_base_key, digits = 3, caption = "Reliability-corrected SEM: key structural paths and effects")
} else {
  cat("No reliability-corrected parameter table is available because the model did not converge.\n")
}
```

## 13. Local-misfit diagnostics from the base reliability-corrected model

```{r rc-base-modification-indices}
rc_base_mi_samefactor <- tibble::tibble(
  lhs = character(),
  op = character(),
  rhs = character(),
  mi = numeric(),
  epc = numeric(),
  sepc.all = numeric()
)

if (is_converged(fit_sem_rc_base)) {
  rc_base_mi_samefactor <- lavaan::modindices(fit_sem_rc_base, sort. = TRUE) %>%
    tibble::as_tibble() %>%
    dplyr::filter(
      op == "~~",
      lhs != rhs,
      mi >= 10,
      (
        (lhs %in% wellbeing_items & rhs %in% wellbeing_items) |
          (lhs %in% belonging_items & rhs %in% belonging_items)
      )
    ) %>%
    dplyr::select(lhs, op, rhs, mi, epc, sepc.all) %>%
    dplyr::slice_head(n = 15)
}

if (nrow(rc_base_mi_samefactor) > 0) {
  knitr::kable(rc_base_mi_samefactor, digits = 3, caption = "Base RC model: top within-construct residual-covariance modification indices")
} else {
  cat("No large same-construct residual-covariance modification indices were found at MI >= 10.\n")
}
```

## 14. Refined reliability-corrected SEM with `Q24a ~~ Q24b`

```{r rc-q24ab-model}
sem_rc_q24ab_lines <- c(
  "family_lat =~ 1*family_support",
  "peer_lat =~ 1*peer_support",
  paste0("family_support ~~ ", format(err_var_family, digits = 8, nsmall = 6), "*family_support"),
  paste0("peer_support ~~ ", format(err_var_peer, digits = 8, nsmall = 6), "*peer_support"),
  "",
  "wellbeing =~ Q24a + Q24b + Q24c + Q24d",
  "school_belonging =~ Q34a + Q34b + Q34c + Q34d + Q34e",
  "",
  "Q24a ~~ Q24b",
  "",
  "school_belonging ~ a1*family_lat + a2*peer_lat",
  "wellbeing ~ b*school_belonging + c1*family_lat + c2*peer_lat",
  "",
  "family_lat ~~ peer_lat",
  "",
  "ind_family := a1*b",
  "ind_peer   := a2*b",
  "total_family := c1 + (a1*b)",
  "total_peer   := c2 + (a2*b)"
)

sem_rc_q24ab_model <- paste(sem_rc_q24ab_lines, collapse = "\n")

fit_sem_rc_q24ab <- tryCatch(
  lavaan::sem(
    model = sem_rc_q24ab_model,
    data = data_sem,
    estimator = lavaan_estimator,
    fixed.x = FALSE,
    missing = missing_handling
  ),
  error = function(e) e
)
```

```{r rc-q24ab-fit}
sem_rc_q24ab_fit <- extract_fit_table(fit_sem_rc_q24ab)
knitr::kable(sem_rc_q24ab_fit, digits = 4, caption = "Refined reliability-corrected SEM fit indices")
```

```{r rc-q24ab-key}
sem_rc_q24ab_key <- extract_key_paths_rc(fit_sem_rc_q24ab)

if (nrow(sem_rc_q24ab_key) > 0) {
  knitr::kable(sem_rc_q24ab_key, digits = 3, caption = "Refined reliability-corrected SEM: key structural paths and effects")
} else {
  cat("No refined-model parameter table is available because the model did not converge.\n")
}
```

## 15. Corrected comparison: base RC model versus refined RC model

```{r rc-model-summary}
rc_model_summary <- dplyr::bind_rows(
  extract_model_row(fit_sem_rc_base, "RC base model"),
  extract_model_row(fit_sem_rc_q24ab, "RC refined model (Q24a ~~ Q24b)")
)

knitr::kable(rc_model_summary, digits = 4, caption = "Model-wise fit summary: base and refined reliability-corrected models")
```

```{r rc-diff-test}
rc_diff_test <- extract_nested_test(
  fit_restricted = fit_sem_rc_base,
  fit_full = fit_sem_rc_q24ab,
  restricted_label = "RC base model",
  full_label = "RC refined model (Q24a ~~ Q24b)"
)

knitr::kable(rc_diff_test, digits = 4, caption = "Difference test: base RC model versus refined RC model")
```

## 16. Fit-change table for the RC comparison

```{r fit-delta-table}
base_row <- rc_model_summary %>%
  dplyr::filter(model == "RC base model")

refined_row <- rc_model_summary %>%
  dplyr::filter(model == "RC refined model (Q24a ~~ Q24b)")

fit_delta_table <- tibble::tibble(
  metric = c("cfi_robust", "tli_robust", "rmsea_robust", "srmr", "aic", "bic"),
  base_model = c(
    base_row$cfi_robust,
    base_row$tli_robust,
    base_row$rmsea_robust,
    base_row$srmr,
    base_row$aic,
    base_row$bic
  ),
  refined_model = c(
    refined_row$cfi_robust,
    refined_row$tli_robust,
    refined_row$rmsea_robust,
    refined_row$srmr,
    refined_row$aic,
    refined_row$bic
  )
) %>%
  dplyr::mutate(
    delta_refined_minus_base = refined_model - base_model
  )

knitr::kable(fit_delta_table, digits = 4, caption = "Change in fit indices from the base RC model to the refined RC model")
```

## 17. Path stability check

```{r path-stability}
base_paths_for_compare <- sem_rc_base_key %>%
  dplyr::mutate(path_id = dplyr::if_else(op == ":=", lhs, paste(lhs, op, rhs))) %>%
  dplyr::select(path_id, std.all, pvalue) %>%
  dplyr::rename(
    std_all_base = std.all,
    pvalue_base = pvalue
  )

refined_paths_for_compare <- sem_rc_q24ab_key %>%
  dplyr::mutate(path_id = dplyr::if_else(op == ":=", lhs, paste(lhs, op, rhs))) %>%
  dplyr::select(path_id, std.all, pvalue) %>%
  dplyr::rename(
    std_all_refined = std.all,
    pvalue_refined = pvalue
  )

path_stability_table <- dplyr::full_join(
  base_paths_for_compare,
  refined_paths_for_compare,
  by = "path_id"
) %>%
  dplyr::mutate(
    delta_std_all = std_all_refined - std_all_base
  )

knitr::kable(path_stability_table, digits = 4, caption = "Path stability: standardized effects before and after adding Q24a ~~ Q24b")
```

## 18. Residual-covariance diagnostics after the refinement

```{r rc-q24ab-modification-indices}
rc_q24ab_mi_samefactor <- tibble::tibble(
  lhs = character(),
  op = character(),
  rhs = character(),
  mi = numeric(),
  epc = numeric(),
  sepc.all = numeric()
)

if (is_converged(fit_sem_rc_q24ab)) {
  rc_q24ab_mi_samefactor <- lavaan::modindices(fit_sem_rc_q24ab, sort. = TRUE) %>%
    tibble::as_tibble() %>%
    dplyr::filter(
      op == "~~",
      lhs != rhs,
      mi >= 10,
      (
        (lhs %in% wellbeing_items & rhs %in% wellbeing_items) |
          (lhs %in% belonging_items & rhs %in% belonging_items)
      )
    ) %>%
    dplyr::select(lhs, op, rhs, mi, epc, sepc.all) %>%
    dplyr::slice_head(n = 15)
}

if (nrow(rc_q24ab_mi_samefactor) > 0) {
  knitr::kable(rc_q24ab_mi_samefactor, digits = 3, caption = "Top within-construct residual-covariance modification indices after adding Q24a ~~ Q24b")
} else {
  cat("No large same-construct residual-covariance modification indices remained at MI >= 10 after adding Q24a ~~ Q24b.\n")
}
```

## 19. SEM path diagrams

This section creates two figures from the **refined reliability-corrected SEM**:

1. a **journal-style structural diagram**  
2. a **full measurement-model diagram** 

```{r sem-diagram-files}
diagram_dir <- file.path("outputs_final_corrected", "figures")
if (!dir.exists(diagram_dir)) dir.create(diagram_dir, recursive = TRUE)

diagram_manifest <- dplyr::bind_rows(
  save_structural_path_diagram(
    fit = fit_sem_rc_q24ab,
    file_stub = "Figure1_refined_sem_structural",
    diagram_dir = diagram_dir,
    width_px = 3800,
    height_px = 2200,
    width_in = 12.5,
    height_in = 7.2
  ),
  save_full_measurement_diagram(
    fit = fit_sem_rc_q24ab,
    file_stub = "FigureS1_refined_sem_full_measurement",
    diagram_dir = diagram_dir,
    width_px = 4200,
    height_px = 3200,
    width_in = 14,
    height_in = 10
  )
)

knitr::kable(diagram_manifest, caption = "Diagram files created from the refined SEM")
```

```{r show-structural-diagram, echo=FALSE, fig.cap="Journal-style structural SEM diagram for the refined reliability-corrected model. Latent constructs are displayed as ellipses, standardized coefficients are shown on all structural paths, and R² values are displayed for the endogenous constructs. The residual covariance Q24a ~~ Q24b is estimated in the fitted model but omitted from the figure for visual clarity."}
struct_png <- file.path(diagram_dir, "Figure1_refined_sem_structural.png")

if (file.exists(struct_png)) {
  knitr::include_graphics(struct_png)
} else {
  cat("The structural SEM diagram was not created.\n")
}
```

```{r show-full-measurement-diagram, echo=FALSE, fig.cap="Full measurement-model SEM diagram for the refined reliability-corrected model."}
full_png <- file.path(diagram_dir, "FigureS1_refined_sem_full_measurement.png")

if (file.exists(full_png)) {
  knitr::include_graphics(full_png)
} else {
  cat("The full measurement-model SEM diagram was not created.\n")
}
```

## 20. Recommendation logic

```{r recommendation-logic}
base_cfi   <- base_row$cfi_robust
base_tli   <- base_row$tli_robust
base_rmsea <- base_row$rmsea_robust
base_srmr  <- base_row$srmr
base_aic   <- base_row$aic
base_bic   <- base_row$bic

ref_cfi   <- refined_row$cfi_robust
ref_tli   <- refined_row$tli_robust
ref_rmsea <- refined_row$rmsea_robust
ref_srmr  <- refined_row$srmr
ref_aic   <- refined_row$aic
ref_bic   <- refined_row$bic

lrt_p <- rc_diff_test$p_value[1]

recommendation_table <- tibble::tibble(
  criterion = c(
    "Refined model converged",
    "Robust CFI improved",
    "Robust TLI improved",
    "Robust RMSEA decreased",
    "SRMR decreased",
    "AIC decreased",
    "BIC decreased",
    "Difference test favours refined model (p < .05)"
  ),
  result = c(
    is_converged(fit_sem_rc_q24ab),
    !is.na(ref_cfi)   && !is.na(base_cfi)   && ref_cfi > base_cfi,
    !is.na(ref_tli)   && !is.na(base_tli)   && ref_tli > base_tli,
    !is.na(ref_rmsea) && !is.na(base_rmsea) && ref_rmsea < base_rmsea,
    !is.na(ref_srmr)  && !is.na(base_srmr)  && ref_srmr < base_srmr,
    !is.na(ref_aic)   && !is.na(base_aic)   && ref_aic < base_aic,
    !is.na(ref_bic)   && !is.na(base_bic)   && ref_bic < base_bic,
    !is.na(lrt_p) && lrt_p < 0.05
  )
)

knitr::kable(recommendation_table, caption = "Rule-based recommendation diagnostics for the refined model")
```

## 21. Model status summary

```{r model-status}
model_status <- tibble::tibble(
  model = c(
    "CFA",
    "Baseline partial mediation",
    "Full mediation",
    "RC base model",
    "RC refined model (Q24a ~~ Q24b)"
  ),
  converged = c(
    is_converged(fit_cfa),
    is_converged(fit_sem_partial),
    is_converged(fit_sem_full),
    is_converged(fit_sem_rc_base),
    is_converged(fit_sem_rc_q24ab)
  ),
  note = c(
    model_note(fit_cfa),
    model_note(fit_sem_partial),
    model_note(fit_sem_full),
    model_note(fit_sem_rc_base),
    model_note(fit_sem_rc_q24ab)
  )
)

knitr::kable(model_status, caption = "Model convergence summary")
```

## 22. Export key outputs

```{r export-results}
output_dir <- "outputs_final_corrected"
if (!dir.exists(output_dir)) dir.create(output_dir)

write.csv(item_dictionary, file.path(output_dir, "item_dictionary.csv"), row.names = FALSE)
write.csv(descriptives, file.path(output_dir, "descriptive_statistics.csv"), row.names = FALSE)
write.csv(missingness, file.path(output_dir, "missingness_table.csv"), row.names = FALSE)
write.csv(reliability_table, file.path(output_dir, "scale_reliability.csv"), row.names = FALSE)
write.csv(cfa_fit_table, file.path(output_dir, "cfa_fit_indices.csv"), row.names = FALSE)
write.csv(sem_partial_fit, file.path(output_dir, "sem_partial_fit_indices.csv"), row.names = FALSE)
write.csv(sem_full_fit, file.path(output_dir, "sem_full_fit_indices.csv"), row.names = FALSE)
write.csv(sem_rc_base_fit, file.path(output_dir, "sem_rc_base_fit_indices.csv"), row.names = FALSE)
write.csv(sem_rc_q24ab_fit, file.path(output_dir, "sem_rc_q24ab_fit_indices.csv"), row.names = FALSE)
write.csv(reliability_correction_table, file.path(output_dir, "reliability_correction_constants.csv"), row.names = FALSE)

write.csv(full_vs_partial_models, file.path(output_dir, "full_vs_partial_model_summary.csv"), row.names = FALSE)
write.csv(full_vs_partial_test, file.path(output_dir, "full_vs_partial_difference_test.csv"), row.names = FALSE)

write.csv(rc_model_summary, file.path(output_dir, "rc_model_summary.csv"), row.names = FALSE)
write.csv(rc_diff_test, file.path(output_dir, "rc_difference_test.csv"), row.names = FALSE)
write.csv(fit_delta_table, file.path(output_dir, "rc_fit_delta_table.csv"), row.names = FALSE)
write.csv(path_stability_table, file.path(output_dir, "path_stability_rc_base_vs_q24ab.csv"), row.names = FALSE)
write.csv(recommendation_table, file.path(output_dir, "recommendation_diagnostics.csv"), row.names = FALSE)
write.csv(model_status, file.path(output_dir, "model_status_summary.csv"), row.names = FALSE)
write.csv(diagram_manifest, file.path(output_dir, "diagram_manifest.csv"), row.names = FALSE)

if (nrow(cfa_loadings) > 0) {
  write.csv(cfa_loadings, file.path(output_dir, "cfa_loadings.csv"), row.names = FALSE)
}
if (nrow(sem_partial_key) > 0) {
  write.csv(sem_partial_key, file.path(output_dir, "sem_partial_key_paths.csv"), row.names = FALSE)
}
if (nrow(sem_rc_base_key) > 0) {
  write.csv(sem_rc_base_key, file.path(output_dir, "sem_rc_base_key_paths.csv"), row.names = FALSE)
}
if (nrow(sem_rc_q24ab_key) > 0) {
  write.csv(sem_rc_q24ab_key, file.path(output_dir, "sem_rc_q24ab_key_paths.csv"), row.names = FALSE)
}
if (nrow(rc_base_mi_samefactor) > 0) {
  write.csv(rc_base_mi_samefactor, file.path(output_dir, "rc_base_modification_indices_samefactor.csv"), row.names = FALSE)
}
if (nrow(rc_q24ab_mi_samefactor) > 0) {
  write.csv(rc_q24ab_mi_samefactor, file.path(output_dir, "rc_q24ab_modification_indices_samefactor.csv"), row.names = FALSE)
}

output_files <- tibble::tibble(
  exported_file = list.files(output_dir, full.names = TRUE, recursive = TRUE)
)

knitr::kable(output_files, caption = "Files written to the outputs_final_corrected folder")
```

## 24. Session information

```{r session-info}
sessionInfo()
```